#!/usr/bin/env nextflow

params.input
params.output_dir = 'rnavar_out'
params.fasta
params.dict
params.known_site1
params.known_site2
params.tbi1
params.tbi2
params.interval_bed
params.bam_intervals
params.vcf_intervals
params.n_xargs = 4

log.info """\
    BILL RNAVAR   P I P E L I N E
    ===================================
    input     : ${params.input}
    fasta     : ${params.fasta}
    dict      : ${params.dict}
    """
    .stripIndent()

process removeDuplicate {
  conda '/home/gjsx/micromamba/envs/gatk'
  cpus 16
  memory 36.GB
  tag "${basename}"
  
  input:
    tuple val(basename), path(input_bam)
  
  output:
    val basename
    path '*.dedup.bam'

  """
  samtools index -@16 $input_bam
  
  gatk MarkDuplicates --java-options "-Xmx36g" \
    --INPUT ${input_bam} \
    --OUTPUT ${basename}.dedup.bam \
    --METRICS_FILE dedup.metrics --TMP_DIR . -REMOVE_DUPLICATES true \
    --ASSUME_SORTED true --VALIDATION_STRINGENCY LENIENT
  """
}

process splitNcigar {
  conda '/home/gjsx/micromamba/envs/gatk'
  cpus 16
  memory 144.GB
  tag "${basename}"
  
  input:
  val basename
  path input_bam
  path bam_intervals
  
  output:
  val basename
  path '*.splited.bam'
  
  """
  samtools index -@16 $input_bam
  
  xargs -P 4 -a ${bam_intervals} \
  -n 2 bash myspliter.sh ${input_bam}
  
  samtools merge -o ${basename}.splited.bam *scattered.bam
  """
}

process recalibrateBQSR {
  conda '/home/gjsx/micromamba/envs/gatk'
  cpus 16
  memory 12.GB
  tag "${basename}"
  
  input:
  val basename
  path input_bam
  path interval_bed
  path fasta
  path dict
  path known_site1
  path known_site2
  path tbi1
  path tbi2
  
  output:
  path '*recal.bam'
  
  """
  samtools index -@16 ${input_bam}
  
  samtools faidx ${fasta}
  
  gatk --java-options "-Xmx12g" BaseRecalibrator  \
        --input ${input_bam} \
        --output recal.table \
        --reference ${fasta} \
        --intervals ${interval_bed} \
        --known-sites ${known_site1} --known-sites ${known_site2} \
        --tmp-dir . \
        --use-original-qualities
  
  gatk --java-options "-Xmx12g" ApplyBQSR \
        --input $input_bam \
        --output ${basename}.recal.bam \
        --reference ${fasta} \
        --bqsr-recal-file recal.table \
        --intervals ${interval_bed} \
        --tmp-dir . \
        --use-original-qualities --add-output-sam-program-record
  """
}

process haploCall {
  conda '/home/gjsx/micromamba/envs/gatk'
  cpus 40
  memory 144.GB
  
  input:
  path input_bam
  path fasta
  path dict
  path vcf_intervals
  
  output:
  path 'merged.bam*'
  path 'haploCall.vcf.gz*'
  path 'filtered.vcf.gz*'
  
  publishDir "${params.output_dir}", mode: 'move'
  
  """
  samtools faidx ${fasta}
  
  samtools merge -@16 --write-index \
  --reference ${fasta} \
  -o merged.bam ${input_bam}
  
  xargs -P 4 -a ${vcf_intervals} \
  mycaller.sh merged.bam
  
  ls -1 *.vcf.gz > vcf_list
  
  gatk --java-options "-Xmx36g" MergeVcfs \
        --INPUT vcf_list \
        --OUTPUT haploCall.vcf.gz \
        --SEQUENCE_DICTIONARY ${dict} \
        --TMP_DIR .
  
  gatk --java-options "-Xmx36G" VariantFiltration \
        --variant haploCall.vcf.gz \
        --output filtered.vcf.gz \
        --reference ${fasta} \
        --tmp-dir . \
        --window 35 --cluster 3 --filter-name "FS" --filter "FS > 30.0" \
        --filter-name "QD" --filter "QD < 2.0"
  """
}

workflow {
  Channel
    .fromPath(params.input)
    .splitCsv(header: true)
    .set{ ch_input }
  bam_intervals_ch = Channel.fromPath(params.bam_intervals)
  vcf_intervals_ch = Channel.fromPath(params.vcf_intervals)
  fasta_ch = Channel.fromPath(params.fasta)
  dict_ch = Channel.fromPath(params.dict)
  site1_ch = Channel.fromPath(params.known_site1)
  site2_ch = Channel.fromPath(params.known_site2)
  bed_ch = Channel.fromPath(params.interval_bed)
  tbi1_ch = Channel.fromPath(params.tbi1)
  tbi2_ch = Channel.fromPath(params.tbi2)
  
  dedup_ch = removeDuplicate(ch_input)
  splited_ch = splitNcigar(dedup_ch, bam_intervals_ch.first())
  
  recal_ch = recalibrateBQSR(
    splited_ch, bed_ch.first(), fasta_ch.first(), dict_ch.first(),
    site1_ch.first(), site2_ch.first(), tbi1_ch.first(), tbi2_ch.first()
    )
    
  haploCall(recal_ch.collect(), fasta_ch.first(), dict_ch.first(),
  vcf_intervals_ch.first())
}
